library(foreign)
library(haven) # read dta 13
library(sandwich) # robust SE
library(SemiParBIVProbit) # Heckman two-stage
library(lmtest) 
library(epiDisplay)
library(stargazer) # printing tables
library(ggplot2) # plots
library(GGally) # more plots
library(dplyr) # efficiently manipulate dataframe
library(Zelig) # rare events logit
library(pROC)
library(mfx)

# Set working directory to Code Location

load("replication_dataset.RData")

######################################
## Replicating Hall and Ura Table 1 ##
######################################

## GLM will get us the coefficients.  We will then need to get robust standard errors.
tab1.floor <- glm(strike ~ FloorMedSup + time + I(time^2) + I(time^3), data = replication_dataset, family = binomial(link = "logit"))
tab1.floor.se <- sqrt(diag(vcovHC(tab1.floor, type="HC0"))) # SEs

## Repeating the above for senate filibuster and party gate keeper as measures of majority support
tab1.filibuster <- glm(strike ~ FilibusterSup + time + I(time^2) + I(time^3), data=replication_dataset, family=binomial(link="logit"))
tab1.filibuster.se <- sqrt(diag(vcovHC(tab1.filibuster, type="HC0")))

## Party Gate Keeper as critical player
tab1.party <- glm(strike ~ PartyGateSup + time + I(time^2) + I(time^3), data = replication_dataset, family = binomial(link="logit"))
tab1.party.se <- sqrt(diag(vcovHC(tab1.party, type = "HC0")))


## Try replicating all of Table 1; standard errors are not in Robust form
tab1.covar.lab <- c("Pr(Floor Median Support)",
                    "Pr(Filibuster Pivot Support)",
                    "Pr(Party Gatekeeper Support)",
                    "Years without invalidation",
                    "Years without invalidation$^2$", 
                    "Years without invalidation$^3$",
                    "Constant")

stargazer(tab1.floor, tab1.filibuster, tab1.party,
          title = "Testing Judicial Majoritanism",
          out = "tables/tab1_replication.tex",
          label = "tab:1",
          column.labels = c("Floor Median Model", "Senate Filibuster Model", "Party Gatekeeping Model"),
          covariate.labels = tab1.covar.lab,
          dep.var.labels.include = F,
          dep.var.caption = "Logistic Regression for US Supreme court invalidating important federal statutes",
          se = list(tab1.floor.se, tab1.filibuster.se, tab1.party.se),
          digits = 2,
          font.size = "footnotesize",
          star.cutoffs = 0.05,
          notes = paste("N =", nrow(tab1.party$model)),
          notes.append = F,
          keep.stat = c("n", "ll", "wald", "chi2"))



########################################
## Replicating Hall and Ura Figure 1a ## 
########################################

# define "one standard deviation above and below the mean"
# plot function takes the variable for measure of majority support, 
# name to be printed on legend, name of PDF, and model to be used
plot.1a <- function(xaxisvar, name, filename, model) {
  
  # setup
  criticalVar <- as.numeric(unlist(replication_dataset[, xaxisvar]))
  lo.value <- mean(criticalVar, na.rm = T) - sd(criticalVar, na.rm = T) 
  hi.value <- mean(criticalVar, na.rm = T) + sd(criticalVar, na.rm = T)
  legend.lab <- c(paste0("Low  (", round(lo.value, 2), ")"),
                  paste0("High (", round(hi.value, 2), ")"))
  
  
  # data frame for plotting
  
  # select variables of interest
  df <- replication_dataset %>%
    dplyr::select(which(colnames(replication_dataset) %in% c("time", xaxisvar))) %>%
    distinct(time)
  
  # set covar to lo value and store predicted probabilities
  df[, 2] <- lo.value
  lo <- predict(model, newdata = df, type  = "response")
  
  # set covar to hi value and store predicted probabilities
  df[, 2] <- hi.value
  hi <- predict(model, newdata = df, type  = "response")
  
  # combine two predicted probabilities
  df <-  df %>%
    arrange(time) %>%
    dplyr::select(time) %>%
    mutate(lo = lo, hi = hi)
  
  # plot
  pdf(paste0("figures/", filename), w = 10)
  plot(df$time, df$lo,
       type = "n", 
       main = "(a) Predicted Probability of Invalidation",
       # main = "Predicted Probability of the Supreme Court invalidating a law\nas the time since enactment or previous invalidation increases\n under different support of Congress",
       ylab = "Predicted Probability of Invalidation",
       xlab = "Years without Invalidation",
       las = 1, 
       bty = "n")
  
  lines(df$time, df$lo, lwd = 3, lty = 2) # low support
  lines(df$time, df$hi, lwd = 3, lty = 1) # high support
  
  legend("topright", legend = legend.lab,
         title = name,
         seg.len = 5,
         lty = c(3, 1), lwd = 3, bty = "n")
  dev.off()
}


plot.1a("FloorMedSup", "Support of Floor Median", "fig1a.pdf", tab1.floor)
plot.1a("FilibusterSup", "Support of Fillbuster Pivot", "fig1a_fillibuster.pdf", tab1.filibuster)
plot.1a("PartyGateSup", "Support of Party Gatekeeper", "fig1a_party.pdf", tab1.party)





########################################
## Replicating Hall and Ura Figure 1b ## 
########################################

# Construct the Model
form.1b.rhs <- c(" ~ time + I(time^2) + I(time^3) + zFloorMedSup") # note we change FloorMed to z

# set values for effect
lo.Maj <- 0
hi.Maj <- 1

z.out.1b <- zelig(as.formula(paste("strike", form.1b.rhs)), model = "logit", data = replication_dataset)

# set up for each t
tticks <- seq(0, 40, 1) #discrete values of time
plot.1b <- matrix(NA, nrow = length(tticks), ncol = 4)
colnames(plot.1b) <- c("pe", "lb", "ub", "time"); rownames(plot.1b) <- tticks

# compute estimate for each t
for(i in seq_along(tticks)){
  x.hi <- setx(z.out.1b, time = tticks[i], zFloorMedSup = hi.Maj) # High Majority Support at time t
  x.lo <- setx(z.out.1b, time = tticks[i], zFloorMedSup = lo.Maj) # low Majority Support
  s.out.1b <- sim(z.out.1b, x = x.lo, x1 = x.hi) # First difference between high and low
  me.1b <- s.out.1b$getqi(qi = "fd", xvalue = "x1") # marginal effects for lo-lo
  plot.1b[i, 1] <-mean(me.1b)
  plot.1b[i, 2:3] <- mean(me.1b) + c(-1, 1)*qnorm(0.975)*sd(me.1b) # CI
  plot.1b[i, 4] <- tticks[i]
}
plot.1b <- data.frame(plot.1b)  # this is what we'll plot

pdf("figures/fig1b.pdf", w = 10)
par(mar = c(5, 6, 3,3))
par(mgp = c(4,1,0))
plot(plot.1b$time, 
     plot.1b$pe, # this goes on y-axis
     ylim = c(-0.020, 0),
     type = "l",
     main = "(b) Marginal Effects of Majority Support",
     xlab = "Years without Invalidation",
     ylab = "Effect of 1 S.D. Increase in Majority Support\n on Predicted Probability of Invalidation",
     las = 1, 
     lwd = 3,
     bty = "n")
lines(plot.1b$time, plot.1b$lb, lwd = 3, lty = 3) # Print Lower CI
lines(plot.1b$time, plot.1b$ub, lwd = 3, lty = 3) # Print Upper CI
abline(h = 0, col = "indianred")
dev.off()



#######################################
## Replicating Hall and Ura Figure 2 ## 
#######################################

# Construct the Model
form.2.rhs <- c(" ~ time + I(time^2) + I(time^3) + FloorMedSup")
z.2 <- zelig(as.formula(paste("challenge", form.2.rhs)), model = "logit", data = replication_dataset)

# set up for each t
floorticks <- seq(0, 1, 0.05) #discrete values of time
plot.2 <- matrix(NA, nrow = length(floorticks), ncol = 4)
colnames(plot.2) <- c("pe", "lb", "ub", "floor"); rownames(plot.2) <- floorticks

# compute estimate for each t
for(i in seq_along(floorticks)){
  x.2 <- setx(z.2, FloorMedSup = floorticks[i], time = 1) # High Majority Support at time t
  s.2 <- Zelig::sim(z.2, x = x.2) # First difference between high and low
  me.2 <- s.2$getqi(qi = "ev") # marginal effects for lo-lo
  plot.2[i, 1] <- mean(me.2)
  plot.2[i, 2:3] <- mean(me.2) + c(-1, 1)*qnorm(0.975)*sd(me.2) # CI why not divide by sqrt(n)?
  plot.2[i, 4] <- floorticks[i]
}
plot.2 <- as.data.frame(plot.2)

# plot
pdf("figures/fig2.pdf", w = 10)
plot(plot.2$floor,
     plot.2$pe, # the DV of interest
     type = "l", 
     main = "",
     xlab = "Majority Support",
     ylab = "Predicted Probability of Challenge",
     ylim = c(0, 0.15),
     xaxt = "n",
     lwd = 3,
     las = 1,
     bty = "n")
axis(1, seq(0, 1, 0.25))
lines(plot.2$floor, plot.2$lb, lwd = 3, lty = 3) # Print Lower CI
lines(plot.2$floor, plot.2$ub, lwd = 3, lty = 3) # Print Upper CI
dev.off()



######################################
## Replicating Hall and Ura Table 2 ##
######################################

# covars same across all models
covars <- paste("time", "time2", "time3", sep = " + ")

## Floor Median Model
tab2.floor.1st <- as.formula(paste("challenge ~", "FloorMedSup +",  covars))
tab2.floor.2nd <- as.formula(paste("strike ~", "FloorMedSup"))
tab2.floor <- SemiParBIVProbit(list(tab2.floor.1st, tab2.floor.2nd), data=replication_dataset, Model="BSS")
summary(tab2.floor) #Replicates coefficients. Stage 2 standard errors off by 0.01


tab2.floor.stage1 <- glm(tab2.floor.1st,data=replication_dataset,family=binomial(link="probit"))
tab2.floor.se <- rep(NA,7)
tab2.floor.se[1:5] <- sqrt(diag(vcovHC(tab2.floor.stage1,type="HC1"))) #replicates Stage 1 standard errors
tab2.floor.se[6:7] <- sqrt(diag(tab2.floor$Vb[6:7,6:7]))
tab2.floor.se

## Senate Filibuter Model
tab2.filibuster.1st <- as.formula(paste("challenge ~", "FilibusterSup +",  covars))
tab2.filibuster.2nd <- as.formula(paste("strike ~", "FilibusterSup"))
tab2.filibuster <- SemiParBIVProbit(list(tab2.filibuster.1st,tab2.filibuster.2nd),data=replication_dataset,Model="BSS")
summary(tab2.filibuster) #Stage 2 SEs off by .01

tab2.filibuster.stage1 <- glm(tab2.filibuster.1st,data=replication_dataset,family=binomial(link="probit"))
tab2.filibuster.se <- rep(NA,7)
tab2.filibuster.se[1:5] <- sqrt(diag(vcovHC(tab2.filibuster.stage1,type="HC1"))) #off by 0.01
tab2.filibuster.se[6:7] <- sqrt(diag(tab2.filibuster$Vb[6:7,6:7]))
tab2.filibuster.se

## Party Gate Keeper Model
tab2.party.1st <- as.formula(paste("challenge ~", "PartyGateSup +",  covars))
tab2.party.2nd <- as.formula(paste("strike ~", "PartyGateSup"))
tab2.party <- SemiParBIVProbit(list(tab2.party.1st,tab2.party.2nd),data=replication_dataset,Model="BSS")
summary(tab2.party) #Stage 2 SEs off by .01

tab2.party.stage1 <- glm(tab2.party.1st,data=replication_dataset,family=binomial(link="probit"))
tab2.party.stage1 <- glm(tab2.filibuster.1st,data=replication_dataset,family=binomial(link="probit"))
tab2.party.se <- rep(NA,7)
tab2.party.se[1:5] <- sqrt(diag(vcovHC(tab2.party.stage1,type="HC1"))) #Correct
tab2.party.se[6:7] <- sqrt(diag(tab2.party$Vb[6:7,6:7]))
tab2.party.se

c(tab2.floor$logLik,tab2.filibuster$logLik,tab2.party$logLik)

## We hand-coded the table as a tex file (see attached file)


#######################################
## Replicating Hall and Ura Figure 3 ##
#######################################


# set values for effect
lo.Maj <- -1
mean.Maj <- 0 # zFloorMedSup is standardized such that SD is 1
hi.Maj <- 1 
lo.agree <- mean(replication_dataset$agree, na.rm = T) - sd(replication_dataset$agree, na.rm = T) 
hi.agree <- mean(replication_dataset$agree, na.rm = T) + sd(replication_dataset$agree, na.rm = T) 
lo.constraint<- mean(replication_dataset$b_con, na.rm = T) - sd(replication_dataset$b_con, na.rm = T) 
hi.constraint <- mean(replication_dataset$b_con, na.rm = T) + sd(replication_dataset$b_con, na.rm = T) 

quadrants <- c("(a) Low Shared Preferences,\n Low Ideological Constraint",
               "(b) High Shared Preferences,\n Low Ideological Constraint",
               "(c) Low Shared Preferences,\n High Ideological Constraint",
               "(d) High Shared Preferences,\n High Ideological Constraint")

# Construct the Model
form.3.rhs <- c(" ~ time + I(time^2) + I(time^3) + agree*b_con*zFloorMedSup") # note we change FloorMed to z
z.out.3 <- Zelig::zelig(as.formula(paste("challenge", form.3.rhs)), model = "logit", data = replication_dataset)


# Lo Lo
x.hi <- setx(z.out.3, agree = lo.agree, b_con = lo.constraint, time = 1, zFloorMedSup = hi.Maj) # predicted value at all values fixed. 
x.lo <- setx(z.out.3, agree = lo.agree, b_con = lo.constraint, time = 1, zFloorMedSup = mean.Maj) # predicted value at all values fixed. 
s.out.3 <- Zelig::sim(z.out.3, x = x.lo, x1 = x.hi)
me.3a <- s.out.3$getqi(qi = "fd", xvalue = "x1") # marginal effects for lo-lo
# Hi Lo
x.hi <- setx(z.out.3, agree = hi.agree, b_con = lo.constraint, time = 1, zFloorMedSup = hi.Maj)
x.lo  <- setx(z.out.3, agree = hi.agree, b_con = lo.constraint, time = 1, zFloorMedSup = mean.Maj)
s.out.3 <-  Zelig::sim(z.out.3, x = x.lo, x1 = x.hi)
me.3b <- s.out.3$getqi(qi = "fd", xvalue = "x1") # marginal effects
# Lo Lo
x.hi <- setx(z.out.3, agree = lo.agree, b_con = hi.constraint, time = 1, zFloorMedSup = hi.Maj)
x.lo  <- setx(z.out.3, agree = lo.agree, b_con = hi.constraint, time = 1, zFloorMedSup = mean.Maj)
s.out.3 <-  Zelig::sim(z.out.3, x = x.lo, x1 = x.hi)
me.3c <- s.out.3$getqi(qi = "fd", xvalue = "x1") # marginal effects
# Hi Lo
x.hi <- setx(z.out.3, agree = hi.agree, b_con = hi.constraint, time = 1, zFloorMedSup = hi.Maj)
x.lo  <- setx(z.out.3, agree = hi.agree, b_con = hi.constraint, time = 1, zFloorMedSup = mean.Maj)
s.out.3 <-  Zelig::sim(z.out.3, x = x.lo, x1 = x.hi)
me.3d <- s.out.3$getqi(qi = "fd", xvalue = "x1") # marginal effects

me.3 <- data.frame(a = me.3a, b = me.3b, c = me.3c, d = me.3d)
me.3.pe <- colMeans(me.3) # marginal effect point estimates
me.3.lb <- me.3.pe - qnorm(0.975)*apply(me.3, 2, sd) # 95 % CI
me.3.ub <- me.3.pe + qnorm(0.975)*apply(me.3, 2, sd)

pdf("figures/fig3.pdf", w = 10)
par(mar = c(3, 6, 3,3))
plot(x = 1:4,
     y = me.3.pe,
     ylim = c(-.04, .04),
     pch = 19,
     las = 1,
     cex = 2,
     xlab = "",
     ylab = "Effect of 1 SD Increase in Majority Support\n on Predicted Probability of Invalidation",
     xaxt = "n",
     bty = "n")
for(i in 1:4) {lines(rep(i, 2), c(me.3.lb[i], me.3.ub[i]), lwd = 3)}
abline(h = 0, lty = 3)
axis(1, 1:4, quadrants, tick = F, cex.axis = 0.8)
dev.off()

########################################
## Replicating Hall and Ura Figure 4 ###
########################################
lo.agree <- mean(replication_dataset$agree, na.rm = T) - sd(replication_dataset$agree, na.rm = T) 
hi.agree <- mean(replication_dataset$agree, na.rm = T) + sd(replication_dataset$agree, na.rm = T) 
lo.constraint<- mean(replication_dataset$b_con, na.rm = T) - sd(replication_dataset$b_con, na.rm = T) 
hi.constraint <- mean(replication_dataset$b_con, na.rm = T) + sd(replication_dataset$b_con, na.rm = T) 

# Construct the Model
form.4.rhs <- c(" ~ time + I(time^2) + I(time^3) + agree*b_con*FloorMedSup")
z.4 <- zelig(as.formula(paste("challenge", form.4.rhs)),model="logit",data=replication_dataset)
summary(z.4)
floorticks <- seq(0, 1, 0.05) #discrete values of FloorMedSup
plot.4 <- matrix(NA, nrow = length(floorticks), ncol = 4)
colnames(plot.4) <- c("pe", "lb", "ub", "floor")

mfx.4 <- function(x) {
  for(i in seq_along(floorticks)){
    x.4 <- setx(z.4, FloorMedSup = floorticks[i], time=1,agree=x[1],b_con=x[2]) #Marginal effects (low shared preferences, low constraint)
    s.4 <-  Zelig::sim(z.4,x=x.4)
    me.4 <- s.4$getqi(qi = "ev", xvalue = "x") # marginal effects
    plot.4[i, 1] <-mean(me.4)
    plot.4[i, 2:3] <- mean(me.4) + c(-1, 1)*qnorm(0.975)*sd(me.4) # CI
    plot.4[i, 4] <- floorticks[i]
  } 
  return(data.frame(plot.4))
}

plot.4a <- mfx.4(c(lo.agree, lo.constraint))
plot.4b <- mfx.4(c(hi.agree, lo.constraint))
plot.4c <- mfx.4(c(lo.agree, hi.constraint))
plot.4d <- mfx.4(c(hi.agree, hi.constraint))
list.4 <- list(plot.4a, plot.4b, plot.4c, plot.4d)


pdf("figures/fig4.pdf", w = 10)
par(mfrow = c(2, 2))
for (i in 1:4){
  plot(list.4[[i]]$floor,
       list.4[[i]]$pe, # the DV of interest
       type = "l", 
       main = quadrants[i],
       xlab = "Majority Support",
       ylab = "Predicted Pr(Court hears challenge to law)",
       ylim = c(0, 0.15),
       xlim = c(.25,1),
       xaxt = "n",
       lwd = 3,
       las = 1,
       bty = "n")
  axis(1, seq(0.25, 1, 0.25))
  lines(list.4[[i]]$floor, list.4[[i]]$lb, lwd = 3, lty = 3) # Print Lower CI
  lines(list.4[[i]]$floor, list.4[[i]]$ub, lwd = 3, lty = 3) # Print Upper CI
}
dev.off()


########################################################
## Alternative formulation of Figure 2: PR Challenge Alt ## 
########################################################

# Construct the Model
form.2.rhs <- c(" ~ time + I(time^2) + I(time^3) + FloorMedSup + I(FloorMedSup^2)")
z.2 <- zelig(as.formula(paste("challenge", form.2.rhs)), model = "logit", data = replication_dataset)

summary(z.2)

# set up for each t
floorticks <- seq(0, 1, 0.05) #discrete values of time
plot.2 <- matrix(NA, nrow = length(floorticks), ncol = 4)
colnames(plot.2) <- c("pe", "lb", "ub", "floor"); rownames(plot.2) <- floorticks

# compute estimate for each t
for(i in seq_along(floorticks)){
  x.2 <- setx(z.2, FloorMedSup = floorticks[i], time = 1) # High Majority Support at time t
  s.2 <- Zelig::sim(z.2, x = x.2) # First difference between high and low
  me.2 <- s.2$getqi(qi = "ev") # marginal effects for lo-lo
  plot.2[i, 1] <- mean(me.2)
  plot.2[i, 2:3] <- mean(me.2) + c(-1, 1)*qnorm(0.975)*sd(me.2) # CI
  plot.2[i, 4] <- floorticks[i]
}
plot.2 <- as.data.frame(plot.2)

# plot
pdf("figures/fig2-alt.pdf", w = 10)
plot(plot.2$floor,
     plot.2$pe, # the DV of interest
     type = "l", 
     main = "",
     xlab = "Majority Support",
     ylab = "Predicted Pr(Court hears challenge to law)",
     ylim = c(0, 0.15),
     xaxt = "n",
     lwd = 3,
     las = 1,
     bty = "n")
axis(1, seq(0, 1, 0.25))
lines(plot.2$floor, plot.2$lb, lwd = 3, lty = 3) # Print Lower CI
lines(plot.2$floor, plot.2$ub, lwd = 3, lty = 3) # Print Upper CI
dev.off()

fig2.main.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + FloorMedSup,family=binomial(link="logit"),data=replication_dataset)
fig2.alt.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + FloorMedSup + I(FloorMedSup^2),family=binomial(link="logit"),data=replication_dataset)
fig2.alt2.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + FloorMedSup + I(FloorMedSup^2) + I(FloorMedSup^3),family=binomial(link="logit"),data=replication_dataset)
summary(fig2.alt.fit)
prob <- predict(fig2.main.fit,type=c("response"))
replication_dataset_subset <- subset(replication_dataset[!is.na(replication_dataset$FloorMedSup),])
replication_dataset_subset$fig2.main.prob <- prob

lrtest(fig2.main.fit,fig2.alt.fit)
lrtest(fig2.main.fit,fig2.alt2.fit)  

# compute ROC
library(pROC)
g <- roc(challenge ~ fig2.main.prob, data = replication_dataset_subset)
prob <- predict(alt.fit,type=c("response"))
replication_dataset_subset$fig2.alt.prob <- prob
f <- roc(challenge ~ fig2.alt.prob, data = replication_dataset_subset)
roc.test(g,f)


########################################################
## Cubic formulation of Figure 2: PR Challenge Alt ## 
########################################################

# Construct the Model
form.2.rhs <- c(" ~ time + I(time^2) + I(time^3) + FloorMedSup + I(FloorMedSup^2) + I(FloorMedSup^3)")
z.2 <- zelig(as.formula(paste("challenge", form.2.rhs)), model = "logit", data = replication_dataset)

summary(z.2)

# set up for each t
floorticks <- seq(0, 1, 0.05) #discrete values of time
plot.2 <- matrix(NA, nrow = length(floorticks), ncol = 4)
colnames(plot.2) <- c("pe", "lb", "ub", "floor"); rownames(plot.2) <- floorticks

# compute estimate for each t
for(i in seq_along(floorticks)){
  x.2 <- setx(z.2, FloorMedSup = floorticks[i], time = 1) # High Majority Support at time t
  s.2 <- Zelig::sim(z.2, x = x.2) # First difference between high and low
  me.2 <- s.2$getqi(qi = "ev") # marginal effects for lo-lo
  plot.2[i, 1] <- mean(me.2)
  plot.2[i, 2:3] <- mean(me.2) + c(-1, 1)*qnorm(0.975)*sd(me.2) # CI
  plot.2[i, 4] <- floorticks[i]
}
plot.2 <- as.data.frame(plot.2)

##########################
## Figure 4 -- Alternative Specification
##########################
lo.agree <- mean(replication_dataset$agree, na.rm = T) - sd(replication_dataset$agree, na.rm = T) 
hi.agree <- mean(replication_dataset$agree, na.rm = T) + sd(replication_dataset$agree, na.rm = T) 
lo.constraint<- mean(replication_dataset$b_con, na.rm = T) - sd(replication_dataset$b_con, na.rm = T) 
hi.constraint <- mean(replication_dataset$b_con, na.rm = T) + sd(replication_dataset$b_con, na.rm = T) 
quadrants <- c("(a) Low Shared Preferences,\n Low Ideological Constraint",
               "(b) High Shared Preferences,\n Low Ideological Constraint",
               "(c) Low Shared Preferences,\n High Ideological Constraint",
               "(d) High Shared Preferences,\n High Ideological Constraint")

# Construct the Model
form.4.rhs.alt <- c(" ~ time + I(time^2) + I(time^3) + agree*b_con*FloorMedSup*I(FloorMedSup^2)")
z.4 <- zelig(as.formula(paste("challenge", form.4.rhs.alt)),model="logit",data=replication_dataset)
summary(z.4)
floorticks <- seq(0, 1, 0.05) #discrete values of FloorMedSup
plot.4 <- matrix(NA, nrow = length(floorticks), ncol = 4)
colnames(plot.4) <- c("pe", "lb", "ub", "floor")

mfx.4 <- function(x) {
  for(i in seq_along(floorticks)){
    x.4 <- setx(z.4, FloorMedSup = floorticks[i], time=1,agree=x[1],b_con=x[2]) #Marginal effects (low shared preferences, low constraint)
    s.4 <-  Zelig::sim(z.4,x=x.4)
    me.4 <- s.4$getqi(qi = "ev", xvalue = "x") # marginal effects
    plot.4[i, 1] <-mean(me.4)
    plot.4[i, 2:3] <- mean(me.4) + c(-1, 1)*qnorm(0.975)*sd(me.4) # CI
    plot.4[i, 4] <- floorticks[i]
  } 
  return(data.frame(plot.4))
}

plot.4a <- mfx.4(c(lo.agree, lo.constraint))
plot.4b <- mfx.4(c(hi.agree, lo.constraint))
plot.4c <- mfx.4(c(lo.agree, hi.constraint))
plot.4d <- mfx.4(c(hi.agree, hi.constraint))
list.4 <- list(plot.4a, plot.4b, plot.4c, plot.4d)


pdf("figures/fig4-alt.pdf", w = 10)
par(mfrow = c(2, 2))
for (i in 1:4){
  plot(list.4[[i]]$floor,
       list.4[[i]]$pe, # the DV of interest
       type = "l", 
       main = quadrants[i],
       xlab = "Majority Support",
       ylab = "Predicted Pr(Court hears challenge to law)",
       ylim = c(0, 0.15),
       xlim = c(.25,1),
       xaxt = "n",
       lwd = 3,
       las = 1,
       bty = "n")
  axis(1, seq(0.25, 1, 0.25))
  lines(list.4[[i]]$floor, list.4[[i]]$lb, lwd = 3, lty = 3) # Print Lower CI
  lines(list.4[[i]]$floor, list.4[[i]]$ub, lwd = 3, lty = 3) # Print Upper CI
}
dev.off()

fig4.main.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*b_con*FloorMedSup,family=binomial(link="logit"),data=replication_dataset)
fig4.cubic.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*b_con*FloorMedSup*I(FloorMedSup^2),family=binomial(link="logit"),data=replication_dataset)
fig4.quadratic.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*b_con*FloorMedSup + agree*b_con*I(FloorMedSup^2),family=binomial(link="logit"),data=replication_dataset)

summary(fig4.cubic.fit)
summary(fig4.quadratic.fit)

#LR test
library(epiDisplay)
lrtest(fig4.main.fit,fig4.cubic.fit) #P-value 0.03
lrtest(fig4.main.fit,fig4.alt2.fit)

#ROC test
prob <- predict(fig4.main.fit,type=c("response"))
replication_dataset_subset <- subset(replication_dataset[!is.na(replication_dataset$FloorMedSup),])
replication_dataset_subset <- subset(replication_dataset_subset[!is.na(replication_dataset_subset$b_con),])
replication_dataset_subset$fig4.main.prob <- prob

# compute ROC
library(pROC)
g <- roc(challenge ~ fig4.main.prob, data = replication_dataset_subset)
prob <- predict(fig4.cubic.fit,type=c("response"))
replication_dataset_subset$fig4.cubic.prob <- prob
f <- roc(challenge ~ fig4.cubic.prob, data = replication_dataset_subset)
roc.test(g,f)

#AIC

AIC(fig4.main.fit)
AIC(fig4.cubic.fit)
AIC(fig4.quadratic.fit)

#BIC

BIC(fig4.main.fit)
BIC(fig4.cubic.fit)
BIC(fig4.quadratic.fit)

#################################################
## Alternative Figure 4--FilibusterSup Measure ##
#################################################

lo.agree <- mean(replication_dataset$agree, na.rm = T) - sd(replication_dataset$agree, na.rm = T) 
hi.agree <- mean(replication_dataset$agree, na.rm = T) + sd(replication_dataset$agree, na.rm = T) 
lo.constraint<- mean(replication_dataset$b_con, na.rm = T) - sd(replication_dataset$b_con, na.rm = T) 
hi.constraint <- mean(replication_dataset$b_con, na.rm = T) + sd(replication_dataset$b_con, na.rm = T) 
quadrants <- c("(a) Low Shared Preferences,\n Low Ideological Constraint",
               "(b) High Shared Preferences,\n Low Ideological Constraint",
               "(c) Low Shared Preferences,\n High Ideological Constraint",
               "(d) High Shared Preferences,\n High Ideological Constraint")

# Construct the Model
form.4.rhs.alt <- c(" ~ time + I(time^2) + I(time^3) + agree*b_con*FilibusterSup*I(FilibusterSup^2)")
z.4 <- zelig(as.formula(paste("challenge", form.4.rhs.alt)),model="logit",data=replication_dataset)
summary(z.4)
floorticks <- seq(0, 1, 0.05) #discrete values of FilibusterSup
plot.4 <- matrix(NA, nrow = length(floorticks), ncol = 4)
colnames(plot.4) <- c("pe", "lb", "ub", "floor")

mfx.4 <- function(x) {
  for(i in seq_along(floorticks)){
    x.4 <- setx(z.4, FilibusterSup = floorticks[i], time=1,agree=x[1],b_con=x[2]) #Marginal effects (low shared preferences, low constraint)
    s.4 <-  Zelig::sim(z.4,x=x.4)
    me.4 <- s.4$getqi(qi = "ev", xvalue = "x") # marginal effects
    plot.4[i, 1] <-mean(me.4)
    plot.4[i, 2:3] <- mean(me.4) + c(-1, 1)*qnorm(0.975)*sd(me.4) # CI
    plot.4[i, 4] <- floorticks[i]
  } 
  return(data.frame(plot.4))
}

plot.4a <- mfx.4(c(lo.agree, lo.constraint))
plot.4b <- mfx.4(c(hi.agree, lo.constraint))
plot.4c <- mfx.4(c(lo.agree, hi.constraint))
plot.4d <- mfx.4(c(hi.agree, hi.constraint))
list.4 <- list(plot.4a, plot.4b, plot.4c, plot.4d)


pdf("figures/fig4-alt-filibuster.pdf", w = 10)
par(mfrow = c(2, 2))
for (i in 1:4){
  plot(list.4[[i]]$floor,
       list.4[[i]]$pe, # the DV of interest
       type = "l", 
       main = quadrants[i],
       xlab = "Majority Support",
       ylab = "Predicted Pr(Court hears challenge to law)",
       ylim = c(0, 0.15),
       xlim = c(.25,1),
       xaxt = "n",
       lwd = 3,
       las = 1,
       bty = "n")
  axis(1, seq(0.25, 1, 0.25))
  lines(list.4[[i]]$floor, list.4[[i]]$lb, lwd = 3, lty = 3) # Print Lower CI
  lines(list.4[[i]]$floor, list.4[[i]]$ub, lwd = 3, lty = 3) # Print Upper CI
}
dev.off()

#LR test
main.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*b_con*FilibusterSup,family=binomial(link="logit"),data=replication_dataset)
alt.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*b_con*FilibusterSup*I(FilibusterSup^2),family=binomial(link="logit"),data=replication_dataset)
lrtest(main.fit,alt.fit) #P-value 0.03

#################################################
## Alternative Figure 4--AvgSup Measure ##
#################################################

lo.agree <- mean(replication_dataset$agree, na.rm = T) - sd(replication_dataset$agree, na.rm = T) 
hi.agree <- mean(replication_dataset$agree, na.rm = T) + sd(replication_dataset$agree, na.rm = T) 
lo.constraint<- mean(replication_dataset$b_con, na.rm = T) - sd(replication_dataset$b_con, na.rm = T) 
hi.constraint <- mean(replication_dataset$b_con, na.rm = T) + sd(replication_dataset$b_con, na.rm = T) 
quadrants <- c("(a) Low Shared Preferences,\n Low Ideological Constraint",
               "(b) High Shared Preferences,\n Low Ideological Constraint",
               "(c) Low Shared Preferences,\n High Ideological Constraint",
               "(d) High Shared Preferences,\n High Ideological Constraint")
replication_dataset$AvgSup <- (replication_dataset$FloorMedSup + replication_dataset$FilibusterSup + replication_dataset$PartyGateSup)/3
# Construct the Model
form.4.rhs.alt <- c(" ~ time + I(time^2) + I(time^3) + agree*b_con*AvgSup*I(AvgSup^2)")
z.4 <- zelig(as.formula(paste("challenge", form.4.rhs.alt)),model="logit",data=replication_dataset)
summary(z.4)
floorticks <- seq(0, 1, 0.05) #discrete values of AvgSup
plot.4 <- matrix(NA, nrow = length(floorticks), ncol = 4)
colnames(plot.4) <- c("pe", "lb", "ub", "floor")

mfx.4 <- function(x) {
  for(i in seq_along(floorticks)){
    x.4 <- setx(z.4, AvgSup = floorticks[i], time=1,agree=x[1],b_con=x[2]) #Marginal effects (low shared preferences, low constraint)
    s.4 <-  Zelig::sim(z.4,x=x.4)
    me.4 <- s.4$getqi(qi = "ev", xvalue = "x") # marginal effects
    plot.4[i, 1] <-mean(me.4)
    plot.4[i, 2:3] <- mean(me.4) + c(-1, 1)*qnorm(0.975)*sd(me.4) #CI
    plot.4[i, 4] <- floorticks[i]
  } 
  return(data.frame(plot.4))
}

plot.4a <- mfx.4(c(lo.agree, lo.constraint))
plot.4b <- mfx.4(c(hi.agree, lo.constraint))
plot.4c <- mfx.4(c(lo.agree, hi.constraint))
plot.4d <- mfx.4(c(hi.agree, hi.constraint))
list.4 <- list(plot.4a, plot.4b, plot.4c, plot.4d)


pdf("figures/fig4-alt-avg.pdf", w = 10)
par(mfrow = c(2, 2))
for (i in 1:4){
  plot(list.4[[i]]$floor,
       list.4[[i]]$pe, # the DV of interest
       type = "l", 
       main = quadrants[i],
       xlab = "Majority Support",
       ylab = "Predicted Pr(Court hears challenge to law)",
       ylim = c(0, 0.15),
       xlim = c(.25,1),
       xaxt = "n",
       lwd = 3,
       las = 1,
       bty = "n")
  axis(1, seq(0.25, 1, 0.25))
  lines(list.4[[i]]$floor, list.4[[i]]$lb, lwd = 3, lty = 3) # Print Lower CI
  lines(list.4[[i]]$floor, list.4[[i]]$ub, lwd = 3, lty = 3) # Print Upper CI
}
dev.off()

#LR test
avg.main.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*b_con*AvgSup,family=binomial(link="logit"),data=replication_dataset)
avg.cubic.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*b_con*AvgSup*I(AvgSup^2),family=binomial(link="logit"),data=replication_dataset)
avg.quadratic.fit <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*b_con*AvgSup + agree*b_con*I(AvgSup^2),family=binomial(link="logit"),data=replication_dataset)
summary(avg.cubic.fit)
lrtest(avg.main.fit,avg.cubic.fit) #P-value 0.023
lrtest(avg.main.fit,avg.quadratic.fit)
lrtest(avg.quadratic.fit,avg.cubic.fit)

avg.cubic.fit.bailey <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*min_dist*AvgSup*I(AvgSup^2),family=binomial(link="logit"),data=replication_dataset)
summary(avg.cubic.fit.bailey)

avg.cubic.fit.bailey2 <- glm(challenge ~ time + I(time^2) + I(time^3) + agree*hs_cons*AvgSup*I(AvgSup^2),family=binomial(link="logit"),data=replication_dataset)

AIC(avg.cubic.fit.bailey2)
AIC(avg.cubic.fit.bailey)
AIC(avg.cubic.fit)

#new plot: AvgSup, cubic, min_dist
lo.constraint.new <- mean(replication_dataset$min_dist, na.rm = T) - sd(replication_dataset$min_dist, na.rm = T) 
hi.constraint.new <- mean(replication_dataset$min_dist, na.rm = T) + sd(replication_dataset$min_dist, na.rm = T) 

form.4.rhs.alt <- c(" ~ time + I(time^2) + I(time^3) + agree*min_dist*AvgSup*I(AvgSup^2)")
z.4 <- zelig(as.formula(paste("challenge", form.4.rhs.alt)),model="logit",data=replication_dataset)
summary(z.4)
floorticks <- seq(0, 1, 0.05) #discrete values of FloorMedSup
plot.4 <- matrix(NA, nrow = length(floorticks), ncol = 4)
colnames(plot.4) <- c("pe", "lb", "ub", "floor")

mfx.4.new <- function(x) {
  for(i in seq_along(floorticks)){
    x.4 <- setx(z.4, AvgSup = floorticks[i], time=1,agree=x[1],min_dist=x[2]) #Marginal effects (low shared preferences, low constraint)
    s.4 <-  Zelig::sim(z.4,x=x.4)
    me.4 <- s.4$getqi(qi = "ev", xvalue = "x") # marginal effects
    plot.4[i, 1] <-mean(me.4)
    plot.4[i, 2:3] <- mean(me.4) + c(-1, 1)*qnorm(0.975)*sd(me.4) # CI
    plot.4[i, 4] <- floorticks[i]
  } 
  return(data.frame(plot.4))
}

plot.4a <- mfx.4.new(c(lo.agree, lo.constraint.new))
plot.4b <- mfx.4.new(c(hi.agree, lo.constraint.new))
plot.4c <- mfx.4.new(c(lo.agree, hi.constraint.new))
plot.4d <- mfx.4.new(c(hi.agree, hi.constraint.new))
list.4 <- list(plot.4a, plot.4b, plot.4c, plot.4d)


pdf("figures/fig4-alt-avg-mindist.pdf", w = 10)
par(mfrow = c(2, 2))
for (i in 1:4){
  plot(list.4[[i]]$floor,
       list.4[[i]]$pe, # the DV of interest
       type = "l", 
       main = quadrants[i],
       xlab = "Majority Support",
       ylab = "Predicted Pr(Court hears challenge to law)",
       ylim = c(0, 0.15),
       xlim = c(.25,1),
       xaxt = "n",
       lwd = 3,
       las = 1,
       bty = "n")
  axis(1, seq(0.25, 1, 0.25))
  lines(list.4[[i]]$floor, list.4[[i]]$lb, lwd = 3, lty = 3) # Print Lower CI
  lines(list.4[[i]]$floor, list.4[[i]]$ub, lwd = 3, lty = 3) # Print Upper CI
}
dev.off()


##
## Fixed effects
##

sillymodel <- glm(challenge ~ time + I(time^2) + I(time^3) + min_dist*agree*AvgSup*I(AvgSup^2) + as.factor(publiclawnum),data=replication_dataset,family=binomial)
lrtest(avg.cubic.fit.bailey,sillymodel)

logitmfx(challenge ~ time + I(time^2) + I(time^3) + min_dist*agree*AvgSup*I(AvgSup^2) + as.factor(publiclawnum),data=replication_dataset)

#ROC test

# compute ROC
replication_dataset_subset <- replication_dataset
replication_dataset_subset <- replication_dataset_subset[!is.na(replication_dataset_subset$AvgSup),]
replication_dataset_subset <- replication_dataset_subset[!is.na(replication_dataset_subset$min_dist),]
prob <- predict(avg.cubic.fit.bailey,type=c("response"))
replication_dataset_subset$avg.cubic.fit.prob <- prob
g <- roc(challenge ~ avg.cubic.fit.prob, data = replication_dataset_subset)
prob <- predict(sillymodel,type=c("response"))
replication_dataset_subset$sillymodel.prob <- prob
f <- roc(challenge ~ sillymodel.prob, data = replication_dataset_subset)
roc.test(g,f)

# plot ROC
pdf("figures/fig4-alt-roc-fixedeffects.pdf", w = 10)
plot(g,
     lty=2,
     main = "",
     las = 1,
     lwd = 2,
     bty = "n")
lines(f)
dev.off()


